OHI British Columbia | OHI Science | Citation policy
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
library(rgdal)
library(sf)
library(tmap)
source('~/github/ohibc/src/R/common.R') ### an OHIBC specific version of common.R
scenario <- 'v2017'
goal <- 'mar'
dir_git <- path.expand('~/github/ohibc')
dir_goal <- file.path(dir_git, 'prep', goal, scenario)
dir_spatial <- file.path(dir_git, 'prep/_spatial')
dir_data_bc <- file.path(dir_M, 'git-annex/bcprep', '_raw_data')
# library(provRmd); prov_setup()
### set up proj4string options: BC Albers and WGS84
p4s_bcalb <- c('bcalb' = '+init=epsg:3005')
p4s_wgs84 <- c('wgs84' = '+init=epsg:4326')This script prepares layers (employment rates and median income) for Livelihoods goal in British Columbia’s coastal regions.
From Halpern et al. (2014) (OHI California Current):
Livelihood sub-goal: As was done in the global analysis, coastal livelihoods is measured by two equally weighted sub-components, the number of jobs (j), which is a proxy for livelihood quantity, and the median annual household wages (g), which is a proxy for job quality. For jobs and wages we used a no-net loss reference point.
For British Columbia, we do not currently have sector-specific unemployment and wage information. As such we will analyze Livelihoods according to the model:
\(x_{LIV} = \frac{j' + g'}{2}\)
\(j' = \frac{j_c / j_{ref}}{M_c / M_{ref}}\)
where M is each region’s employment rate (1 - unemployment) as a percent at current (c) and reference (ref) time periods, and:
\(g' = \frac{g_c / g_{ref}}{W_c / W_{ref}}\)
where W is each State’s average annual per capita wage at current (c) and reference (ref) time periods.
phi_rasts <- c('tf_rast' = file.path(dir_goal, 'data_explore/t_f_raster_1000m.tif'),
'tb_rast' = file.path(dir_goal, 'data_explore/t_b_raster_1000m.tif'),
'tb_minus_sd_rast' = file.path(dir_goal, 'data_explore/t_b_minus_sd_raster_1000m.tif'))
dir_aq <- file.path(dir_data_bc, 'aquaculture')
reload <- FALSE
if(any(!file.exists(phi_rasts)) | reload) {
ohibc_rast <- raster(file.path(dir_spatial, 'raster/ohibc_rgn_raster_1000m.tif'))
phif_rast <- raster(file.path(dir_aq, 'FishPhiALLConstraints95LT2.tif')) %>%
projectRaster(ohibc_rast)
tf_rast <- exp(7.86 - log(phif_rast) * 5.82)
phib_rast <- raster(file.path(dir_aq, 'BivalvePhiALLConstraints95LT1.tif')) %>%
projectRaster(ohibc_rast)
phib_sd_rast <- raster(file.path(dir_aq, 'Bivalve_spp_Phi_sd.tif')) %>%
projectRaster(ohibc_rast)
# x <- phib_rast - phib_sd_rast
tb_rast <- exp(2.99 - phib_rast * 1.66)
tb_minus_sd_rast <- exp(2.99 - (phib_rast - phib_sd_rast) * 1.66)
writeRaster(tf_rast, phi_rasts['tf_rast'], overwrite = TRUE)
writeRaster(tb_rast, phi_rasts['tb_rast'], overwrite = TRUE)
writeRaster(tb_minus_sd_rast, phi_rasts['tb_minus_sd_rast'], overwrite = TRUE)
}phi_rasts_uncl <- c('tf_rast' = file.path(dir_goal, 'data_explore/t_f_unclipped_1000m.tif'),
'tb_rast' = file.path(dir_goal, 'data_explore/t_b_unclipped_1000m.tif'),
'tb_minus_sd_rast' = file.path(dir_goal, 'data_explore/t_b_minus_sd_unclipped_1000m.tif'))
dir_aq <- file.path(dir_data_bc, 'aquaculture')
reload <- FALSE
if(any(!file.exists(phi_rasts_uncl)) | reload) {
ohibc_rast <- raster(file.path(dir_spatial, 'raster/ohibc_rgn_raster_1000m.tif'))
phif_rast_uncl <- raster(file.path(dir_aq, 'spp_Phi_mean.tif')) %>%
projectRaster(ohibc_rast)
values(phif_rast_uncl)[values(phif_rast_uncl) < 2.0] <- NA
### clip out unreasonably low Phi values
tf_rast_uncl <- exp(7.86 - log(phif_rast_uncl) * 5.82)
phib_rast_uncl <- raster(file.path(dir_aq, 'Bivalve_spp_Phi_mean.tif')) %>%
projectRaster(ohibc_rast)
phib_sd_rast <- raster(file.path(dir_aq, 'Bivalve_spp_Phi_sd.tif')) %>%
projectRaster(ohibc_rast)
tb_rast_uncl <- exp(2.99 - phib_rast_uncl * 1.66)
tb_minus_sd_rast_uncl <- exp(2.99 - (phib_rast_uncl - phib_sd_rast) * 1.66)
writeRaster(tf_rast_uncl, phi_rasts_uncl['tf_rast'], overwrite = TRUE)
writeRaster(tb_rast_uncl, phi_rasts_uncl['tb_rast'], overwrite = TRUE)
writeRaster(tb_minus_sd_rast_uncl, phi_rasts_uncl['tb_minus_sd_rast'], overwrite = TRUE)
}tf_rast <- raster(file.path(dir_goal, 'data_explore/t_f_raster_1000m.tif'))
tb_rast <- raster(file.path(dir_goal, 'data_explore/t_b_raster_1000m.tif'))
tf_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_f_unclipped_1000m.tif'))
tb_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_b_unclipped_1000m.tif'))
ohibc_sf <- sf::read_sf(dir_spatial, 'ohibc_rgn')
bc_land_sf <- sf::read_sf(dir_spatial, 'ohibc_land')
tf_map <- tm_shape(bc_land_sf) +
tm_fill(col = 'grey40', alpha = 1) +
tm_shape(ohibc_sf) +
tm_fill(col = 'grey96', alpha = 1) +
tm_shape(tf_rast_uncl) +
tm_raster(title = 'Growth time\nyrs (35 cm fish)\n(all data)',
palette = 'Greens') +
tm_shape(tf_rast) +
tm_raster(title = 'Growth time\nyears (35 cm fish)',
palette = 'Reds') +
tm_shape(ohibc_sf) +
tm_borders(col = 'grey40', lwd = .25) +
tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))
print(tf_map)
tb_map <- tm_shape(bc_land_sf) +
tm_fill(col = 'grey40', alpha = 1) +
tm_shape(ohibc_sf) +
tm_polygons(col = 'grey96', alpha = 1) +
tm_shape(tb_rast_uncl) +
tm_raster(title = 'Growth time\nyrs (4 cm bivalve)\n(all data)',
palette = 'Greens') +
tm_shape(tb_rast) +
tm_raster(title = 'Growth time\nyears (4 cm bivalve)',
palette = 'Blues') +
tm_shape(ohibc_sf) +
tm_borders(col = 'grey40', lwd = .25) +
tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))
print(tb_map)From the growth time, we calculate biomass production based on the following assumptions and calculations:
tf_rast <- raster(file.path(dir_goal, 'data_explore/t_f_raster_1000m.tif'))
harvest_f_rast <- 2376/tf_rast
tb_rast <- raster(file.path(dir_goal, 'data_explore/t_b_raster_1000m.tif'))
harvest_b_rast <- 130e6/tb_rast
tb_sd_rast <- raster(file.path(dir_goal, 'data_explore/t_b_minus_sd_raster_1000m.tif'))
harvest_b_sd_rast <- 130e6/tb_sd_rast
writeRaster(harvest_f_rast, file.path(dir_goal, 'data_explore/harvest_f_raster_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_rast, file.path(dir_goal, 'data_explore/harvest_b_raster_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_sd_rast, file.path(dir_goal, 'data_explore/harvest_b_sd_raster_1000m.tif'), overwrite = TRUE)
tf_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_f_unclipped_1000m.tif'))
harvest_f_rast_uncl <- 2376/tf_rast_uncl
tb_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_b_unclipped_1000m.tif'))
harvest_b_rast_uncl <- 130e6/tb_rast_uncl
tb_sd_rast_uncl <- raster(file.path(dir_goal, 'data_explore/t_b_minus_sd_unclipped_1000m.tif'))
harvest_b_sd_rast_uncl <- 130e6/tb_sd_rast_uncl
writeRaster(harvest_f_rast_uncl, file.path(dir_goal, 'data_explore/harvest_f_unclipped_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_rast_uncl, file.path(dir_goal, 'data_explore/harvest_b_unclipped_1000m.tif'), overwrite = TRUE)
writeRaster(harvest_b_sd_rast_uncl, file.path(dir_goal, 'data_explore/harvest_b_sd_unclipped_1000m.tif'), overwrite = TRUE)harvest_f_rast <- raster(file.path(dir_goal, 'data_explore/harvest_f_raster_1000m.tif'))
harvest_b_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_raster_1000m.tif'))
harvest_b_sd_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_sd_raster_1000m.tif'))
harvest_f_rast_uncl <- raster(file.path(dir_goal, 'data_explore/harvest_f_unclipped_1000m.tif'), overwrite = TRUE)
harvest_b_rast_uncl <- raster(file.path(dir_goal, 'data_explore/harvest_b_unclipped_1000m.tif'), overwrite = TRUE)
harvest_b_sd_rast_uncl <- raster(file.path(dir_goal, 'data_explore/harvest_b_sd_unclipped_1000m.tif'), overwrite = TRUE)
ohibc_sf <- sf::read_sf(dir_spatial, 'ohibc_rgn')
bc_land_sf <- sf::read_sf(dir_spatial, 'ohibc_land')
harvest_f_map <- tm_shape(bc_land_sf) +
tm_fill(col = 'grey40', alpha = 1) +
tm_shape(ohibc_sf) +
tm_fill(col = 'grey80', alpha = .4) +
tm_shape(harvest_f_rast_uncl) +
tm_raster(title = 'Fish harvest\ntonnes/yr (all pts)',
palette = 'Greens') +
tm_shape(harvest_f_rast) +
tm_raster(title = 'Fish harvest\ntonnes/year') +
tm_shape(ohibc_sf) +
tm_borders(col = 'grey40', lwd = .25) +
tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))
print(harvest_f_map)
harvest_b_map <- tm_shape(bc_land_sf) +
tm_fill(col = 'grey40', alpha = 1) +
tm_shape(ohibc_sf) +
tm_polygons(col = 'grey80', alpha = .4) +
tm_shape(harvest_b_rast_uncl) +
tm_raster(title = 'Bivalve harvest\nunits/yr (all pts)',
palette = 'Greens') +
tm_shape(harvest_b_rast) +
tm_raster(title = 'Bivalve harvest\nunits/yr',
palette = 'Blues') +
tm_shape(ohibc_sf) +
tm_borders(col = 'grey40', lwd = .25) +
tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))
print(harvest_b_map)
harvest_b_sd_map <- tm_shape(bc_land_sf) +
tm_fill(col = 'grey40', alpha = 1) +
tm_shape(ohibc_sf) +
tm_polygons(col = 'grey80', alpha = .4) +
tm_shape(harvest_b_sd_rast_uncl) +
tm_raster(title = 'Bivalve harvest\nunits/yr (-1 sd)',
palette = 'Greens') +
tm_shape(harvest_b_sd_rast) +
tm_raster(title = 'Bivalve harvest\nunits/yr (-1 sd)',
palette = 'Blues') +
tm_shape(ohibc_sf) +
tm_borders(col = 'grey40', lwd = .25) +
tm_legend(bg.alpha = .7, bg.color = 'white', position = c('right', 'top'))
print(harvest_b_sd_map)For each OHIBC region, examine the distribution of production potential for both fish and bivalves. Some possible targets could be based on median production for each region, other quantiles. Due to coarseness of production raster, production hot spots do not line up well with tenure locations, so spatially estimating production potential based on tenures is not likely to be a good method.
harvest_f_rast <- raster(file.path(dir_goal, 'data_explore/harvest_f_raster_1000m.tif'))
### use extract() to get values per region
ohibc_poly <- readOGR(dir_spatial, 'ohibc_rgn')## OGR data source with driver: ESRI Shapefile
## Source: "/home/ohara/github/ohibc/prep/_spatial", layer: "ohibc_rgn"
## with 8 features
## It has 5 fields
## Integer64 fields read as strings: rgn_id
f_values <- raster::extract(harvest_f_rast, ohibc_poly) %>%
lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
setNames(ohibc_poly$rgn_id) %>%
bind_rows(.id = 'rgn_id') %>%
mutate(rgn_id = as.integer(rgn_id),
prod_potential = round(prod_potential, 3))
prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
f_range <- f_values %>%
group_by(rgn_id) %>%
summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)),
min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
unnest() %>%
mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
ungroup()
f_plot <- ggplot(f_values %>%
left_join(get_rgn_names(), by = 'rgn_id')) +
ggtheme_plot() +
geom_histogram(aes(x = prod_potential)) +
facet_wrap(~ rgn_name, scales = 'free_y') +
labs(x = 'Production potential (Finfish tonnes/km^2/year)',
y = 'Number of 1 km^2 cells',
title = 'Finfish')
print(f_plot)write_csv(f_range, file.path(dir_goal, 'data_explore/prod_pot_f_range.csv'))harvest_b_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_raster_1000m.tif'))
b_values <- raster::extract(harvest_b_rast, ohibc_poly) %>%
lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
setNames(ohibc_poly$rgn_id) %>%
bind_rows(.id = 'rgn_id') %>%
mutate(rgn_id = as.integer(rgn_id),
prod_potential = round(prod_potential))
prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
b_range <- b_values %>%
group_by(rgn_id) %>%
summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)),
min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
unnest() %>%
mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
ungroup()
b_plot <- ggplot(b_values %>%
left_join(get_rgn_names(), by = 'rgn_id')) +
ggtheme_plot() +
geom_histogram(aes(x = prod_potential)) +
facet_wrap(~ rgn_name, scales = 'free_y') +
labs(x = 'Production potential (4 cm bivalves/km^2/year)',
y = 'Number of 1 km^2 cells',
title = 'Bivalves')
print(b_plot)write_csv(b_range, file.path(dir_goal, 'data_explore/prod_pot_b_range.csv'))In generating median harvest values for bivalves from the unclipped dataset (to fill the NAs in North Coast, and flesh out a few more datapoints for Strait of Georgia and Haida Gwaii), we want to ensure that the numbers aren’t outrageous. For the regions with data in the clipped dataset, we compare to the unclipped and see that the medians still fall fairly close.
harvest_b_uncl_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_unclipped_1000m.tif'))
b_uncl_values <- raster::extract(harvest_b_uncl_rast, ohibc_poly) %>%
lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
setNames(ohibc_poly$rgn_id) %>%
bind_rows(.id = 'rgn_id') %>%
mutate(rgn_id = as.integer(rgn_id),
prod_potential = round(prod_potential))
prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
b_uncl_range <- b_uncl_values %>%
group_by(rgn_id) %>%
summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)),
min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
unnest() %>%
mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
ungroup()
b_uncl_plot <- ggplot(b_uncl_values %>%
left_join(get_rgn_names(), by = 'rgn_id')) +
ggtheme_plot() +
geom_histogram(aes(x = prod_potential)) +
facet_wrap(~ rgn_name, scales = 'free_y') +
labs(x = 'Production potential (4 cm bivalves/km^2/year) (all data)',
y = 'Number of 1 km^2 cells',
title = 'Bivalves (all data)')
print(b_uncl_plot)write_csv(b_uncl_range, file.path(dir_goal, 'data_explore/prod_pot_b_range_unclipped.csv'))
x <- b_range %>%
select(rgn_id, quants, prob_lvl) %>%
left_join(b_uncl_range %>%
select(rgn_id, quants_unclipped = quants, prob_lvl),
by = c('rgn_id', 'prob_lvl')) %>%
# filter(prob_lvl == .50) %>%
left_join(get_rgn_names())
test_plot <- ggplot(x, aes(x = quants, y = quants_unclipped, label = prob_lvl)) +
geom_abline(intercept = 0, slope = 1, color = 'red') +
geom_hline(data = x %>% filter(is.na(quants) & prob_lvl == 0.5),
aes(yintercept = quants_unclipped, label = rgn_name), alpha = .5) +
geom_point(alpha = .3, show.legend = FALSE) +
geom_point(data = x %>% filter(prob_lvl == 0.5), aes(color = rgn_name)) +
xlim(0, NA) + ylim(0, NA) +
labs(title = 'Median bivalve production, \nall data vs clipped',
x = 'prod. potential clipped',
y = 'prod. potential, all cells')
plotly::ggplotly(test_plot)Since there is such uncertainty around bivalve growth rates let’s examine the impact of subtracting the standard deviation, as a possible lower target. This plot considers unclipped values for all regions.
harvest_b_sd_uncl_rast <- raster(file.path(dir_goal, 'data_explore/harvest_b_sd_unclipped_1000m.tif'))
b_sd_uncl_values <- raster::extract(harvest_b_sd_uncl_rast, ohibc_poly) %>%
lapply(FUN = function(x) data.frame(prod_potential = x)) %>%
setNames(ohibc_poly$rgn_id) %>%
bind_rows(.id = 'rgn_id') %>%
mutate(rgn_id = as.integer(rgn_id),
prod_potential = round(prod_potential))
prob_levels <- c(.01, .02, .05, .10, .25, .50, .75, .90, .95, .98, .99)
b_sd_uncl_range <- b_sd_uncl_values %>%
group_by(rgn_id) %>%
summarize(n_tot = n(), n_phi = sum(!is.na(prod_potential)),
min_prod = min(prod_potential, na.rm = TRUE), max_prod = max(prod_potential, na.rm = TRUE),
quants = list(stats::quantile(prod_potential, probs = prob_levels, na.rm = TRUE))) %>%
unnest() %>%
mutate(prob_lvl = rep(prob_levels, times = 8)) %>%
ungroup()
b_sd_uncl_plot <- ggplot(b_sd_uncl_values %>%
left_join(get_rgn_names(), by = 'rgn_id')) +
ggtheme_plot() +
geom_histogram(aes(x = prod_potential)) +
facet_wrap(~ rgn_name, scales = 'free_y') +
labs(x = 'Production potential (4 cm bivalves/km^2/year) (-1 sd)',
y = 'Number of 1 km^2 cells',
title = "Bivalves (median based on Phi' - 1sd)")
print(b_sd_uncl_plot)write_csv(b_sd_uncl_range, file.path(dir_goal, 'data_explore/prod_pot_b_sd_range_unclipped.csv'))
x <- b_uncl_range %>%
select(rgn_id, quants, prob_lvl) %>%
left_join(b_sd_uncl_range %>%
select(rgn_id, quants_minus_sd = quants, prob_lvl),
by = c('rgn_id', 'prob_lvl')) %>%
# filter(prob_lvl == .50) %>%
left_join(get_rgn_names())
test_plot <- ggplot(x, aes(x = quants, y = quants_minus_sd, label = prob_lvl)) +
geom_abline(intercept = 0, slope = 1, color = 'red') +
geom_point(alpha = .3, show.legend = FALSE) +
geom_point(data = x %>% filter(prob_lvl == 0.5), aes(color = rgn_name)) +
xlim(0, NA) + ylim(0, NA) +
labs(title = 'Bivalve production, -1 sd',
x = 'prod. potential, all cells',
y = "prod. potent, phi' - 1sd")
plotly::ggplotly(test_plot)As aquaculture development targets, we will use MaPP Special Management Zones assigned to aquaculture to represent the goals of the MaPP regions. For non-MaPP regions, we will use currently allocated aquaculture tenures (active/not, developed/not) to indicate development goals.
To determine the area allocated to aquaculture/mariculture for the MaPP regions, we will use SMZ proportions of total region ocean area, taken from MaPP region shapefiles. In some cases (e.g. Haida Gwaii west coast north, table 8.21) it appears some non-SMZ areas are noted “acceptable” for aquaculture. As a conservative estimate, only SMZ areas will be counted.
Note that for North Coast Vancouver Island, the OHIBC region wraps around to the west in addition to the overall MaPP region. In this case, perhaps we want to include existing aquaculture tenure areas? For now, these are excluded.
Do any of these regions allow finfish aquaculture? If so how to allocate the area appropriately…
dir_mapp <- file.path(dir_data_bc, 'mapp/MarineSpatialPlanZones')
mapp_hg <- read_sf(dir_mapp, 'HaidaGwaii_Subregion_Oct_2014')
# "SubRegion" "Type" "Management" "Area_Ha" "Area_Km" "Objective" "Name"
# "Haida_Name" "Shape_Leng" "Shape_Area"
# mapp_hg$Management %>% unique()
# [1] "SMZ - Shellfish Aq" "IUCN - Type Ib" "IUCN - Type IV" "IUCN - Type II"
# [5] "IUCN - Type III" "IUCN - Type V" "IUCN - Type VI" "SMZ - Alt. Energy"
mapp_hg_aq <- mapp_hg %>%
filter(Management == 'SMZ - Shellfish Aq') %>%
mutate(area = st_area(geometry)) %>%
as.data.frame() %>%
dplyr::select(name = Name, area, desc = Management) %>%
mutate(rgn_id = 2)
mapp_nc <- read_sf(dir_mapp, 'NorthCoast_Subregion_June2015')
# "SubRegion" "Type" "Management" "Area_Ha" "Area_Km" "Objective" "AreaDescri"
# "Name" "Unit_No" "Grouping" "ZoneType" "Shape_Leng" "Shape_Area"
# mapp_nc$Management %>% unique()
# [1] "II" "IV"
# [3] "Renewable Energy" "Aquaculture"
# [5] "Tourism and Recreation" "Areas for future planning consideration"
# [7] "Ib" "Cultural"
mapp_nc_aq <- mapp_nc %>%
filter(Management == 'Aquaculture') %>%
mutate(area = st_area(geometry)) %>%
as.data.frame() %>%
dplyr::select(name = Name, area, desc = Management) %>%
mutate(rgn_id = 1)
mapp_cc <- read_sf(dir_mapp, 'CentralCoast_Subregion_Ver9')
# "SubRegion" "Type" "Area_Ha" "Area_Km" "Unit_New" "Unit" "Group_"
# "IUCN_Categ" "ZoneType" "Edited" "Unit_Old" "AreaDecri" "Shape_Leng" "Shape_Area"
# mapp_cc$ZoneType %>% unique()
# [1] NA "Renewable Energy" "Aquaculture"
# [4] "Recreation and Tourism"
mapp_cc_aq <- mapp_cc %>%
filter(ZoneType == 'Aquaculture') %>%
mutate(area = st_area(geometry)) %>%
as.data.frame() %>%
dplyr::select(name = AreaDecri, area, desc = ZoneType) %>%
mutate(rgn_id = 3)
mapp_ncvi <- read_sf(dir_mapp, 'NorthVancouverIsland_Subregion_Oct3_2014')
# "SubRegion" "Type" "Management" "Area_Ha" "Area_Km" "AreaDescri" "Name"
# "Unit_No" "Shape_Leng" "Shape_Area"
# mapp_ncvi$AreaDescri %>% unique()
mapp_ncvi_aq <- mapp_ncvi %>%
filter(str_detect(tolower(AreaDescri), 'aquaculture')) %>%
mutate(area = st_area(geometry)) %>%
as.data.frame() %>%
dplyr::select(name = Name, area, desc = AreaDescri) %>%
mutate(rgn_id = 4)
# [1] "SMZ Cultural/Economic Emphasis Areas are intended to reinforce their high value to First Nations, on a seasonal and year-round basis, for cultural value protection, Aboriginal economic development opportunities, and food security. This value includes con"
# [2] "Significant ecological values due to major upwelling of nutrients creating a rich, diverse marine ecosystem. There are key First Nation cultural/economic interests and local resident scenic values. Safeguarding the integrity of this interaction between"
# [3] "SMZ Recreation/Tourism Emphasis Areas are intended to reinforce their high value to existing commercial tourism operations, particularly during the months of late May to early October. Other uses and activities in SMZ Recreation/Tourism Emphasis Areas sh"
# [4] "SMZ Community Emphasis Areas are intended to reinforce their value for seasonal and year-round uses and activities associated with, required by, and primarily dictated by, adjacent, or nearby communities. The uses and activities in SMZ Community Emphasis"
# [5] "SMZ Shellfish Aquaculture Emphasis Areas are intended to reinforce interest by First Nations in investigation and (if feasible) the development of bottom and off-bottom shellfish aquaculture operations. These areas may be associated with integrated multi"
# [6] "Important species and habitats, including those of cultural importance to First Nations. Significant for seasonal marine harvesting and ecotourism activities by First Nations. It is an important whale and wildlife viewing area. Includes important habitat"
# [7] "Important habitat and species, in particular a significant and unique glass sponge reef formation, which includes a complex ecosystem, enabling a species-rich marine environment that supports the local biodiversity of the area. Contains critical habitat"
# [8] "The area is representative of shallow sill ecosystems of coral fans, sponges. Several deepwater and/or rare species including the gorgonian coral, the soft goblet sponge, the cloud sponge, the townsend eualid shrimp, and the bigmouth sculpin are found at"
# [9] "Considerable cultural modification by First Nations based on use of important marine species and habitat. Ongoing practices and teachings, restoration of First Nations’ cultural resources, and their associated marine species and habitats, and for repa"
# [10] "Important marine species, habitats, First Nations’ cultural resources such as cultural tourism, loxiwe, shell middens, and former seasonal village/resource processing site."
# [11] "High marine recreational values, containing important marine species and habitats including important areas for herring and northern resident killer whales. Area includes First Nations’ cultural resources uses and activities such as cultural tourism, l"
# [12] "A diverse marine ecosystem, with important marine species and habitat. Important recreation and tourism area which includes several scuba diving sites. Includes important areas for humpback and and northern resident killer whales, herring and sea otters."
# [13] "Marine species and habitats including those of cultural importance to First Nations. Connects existing conservation and protection areas and provides network/corridor between the Central Coast and NVI marine plans to assist in conservation and protection"
# [14] "Important marine species and habitats including herring important areas. Protection of representative marine ecosystems at the confluence of three channels supporting rich intertidal species and habitats."
mapp_aq_area_df <- bind_rows(mapp_hg_aq,
mapp_cc_aq,
mapp_nc_aq,
mapp_ncvi_aq) %>%
mutate(area_km2 = round(area / 1e6, 3)) %>%
select(-area)
write_csv(mapp_aq_area_df, file.path(dir_goal, 'data_explore/mapp_aq_smz_areas.csv'))For Strait of Georgia and West Coast Vancouver Island, we will use the total area of designated tenures according to DFO records, split into finfish and shellfish.
Note that for Aristazabal Island OHIBC region, no aquaculture areas (tenures or SMZs) have been identified.
tenure_shps <- list.files(file.path(dir_data_bc, 'dfo_khunter/aquaculture/d2016/TENURES_PNT'),
pattern = '.shp$', full.names = TRUE)
ohibc_rgn <- read_sf(dir_spatial, 'ohibc_rgns_unclipped')
tenures_ohibc <- lapply(tenure_shps, FUN = function(x) {
### x <- tenure_shps[1]
read_sf(dsn = dirname(x), layer = str_replace(basename(x), '.shp$', '')) %>%
mutate(area = st_area(geometry)) %>%
st_intersection(ohibc_rgn) %>%
as.data.frame() %>%
dplyr::select(-geometry)
}) %>%
setNames(str_replace(basename(tenure_shps), '.shp$', '')) %>%
bind_rows(.id = 'filename')
tenure_areas <- tenures_ohibc %>%
setNames(tolower(names(.))) %>%
mutate(spp = ifelse(is.na(sp_lic), sp_, sp_lic),
area_km2 = round(area / 1e6, 3)) %>%
dplyr::select(filename, rgn_id, location, spp, area_km2)
write_csv(tenure_areas, file.path(dir_goal, 'data_explore/dfo_aq_tenure_areas.csv'))DFO data for aquaculture harvests by PFMA for 2011-2015. Harvest values are in tonnes?
aq_dfo_file <- file.path(dir_data_bc, 'dfo_khunter/aquaculture/d2016',
'Aquaculture_production_PFMA 2011-15.xlsx')
data_aq_dfo <- readxl::read_excel(aq_dfo_file, skip = 4) %>%
setNames(tolower(names(.)) %>%
str_replace_all('[^a-z]+', '_'))
### This file has a number of rows for shellfish, then a number of rows for
### finfish. The header row can be detected by looking for non-numeric
### text in the "year" column.
start_fish <- which(data_aq_dfo$year %>% str_detect('[a-zA-Z]'))
aq_shellfish <- data_aq_dfo[1:(min(start_fish) - 1), ] %>%
gather(key = 'species', value = 'harvest', -year, -pfma) %>%
mutate(aq_type = 'shellfish')
finfish_hdr <- as.character(data_aq_dfo[max(start_fish), ]) %>%
tolower() %>%
str_replace_all('[^a-z]+', '_') %>%
.[. != 'na']
aq_finfish <- data_aq_dfo[(max(start_fish) + 1):nrow(data_aq_dfo), 1:length(finfish_hdr)] %>%
setNames(finfish_hdr) %>%
gather(key = 'species', value = 'harvest', -year, -pfma) %>%
mutate(aq_type = 'finfish')
aq_harvest_df <- bind_rows(aq_shellfish, aq_finfish) %>%
filter(!is.na(year))
write_csv(aq_harvest_df, file.path(dir_goal, 'data_explore/dfo_aq_harvest_pfma.csv'))Since the DFO harvest is provided by PFMA, collapse PFMAs to OHIBC regions. We have already done this in the FIS goal, so we will use the same table to assign PFMAs here.
pfma_to_ohibc <- read_csv(file.path(dir_git, 'prep/fis/v2017/int/pfmsa_to_ohibc.csv')) %>%
group_by(pfma_id, rgn_id) %>%
summarize(area_km2 = sum(area_km2)) %>%
group_by(pfma_id) %>%
mutate(prop_area = area_km2 / sum(area_km2)) %>%
select(pfma = pfma_id, rgn_id, prop_area)
aq_harvest_ohibc <- read_csv(file.path(dir_goal, 'data_explore/dfo_aq_harvest_pfma.csv')) %>%
left_join(pfma_to_ohibc, by = 'pfma') %>%
mutate(harvest_rgn = harvest * prop_area) %>%
select(year, species, harvest_rgn, aq_type, rgn_id) %>%
group_by(year, species, aq_type, rgn_id) %>%
summarize(harvest = sum(harvest_rgn)) %>%
ungroup()
write_csv(aq_harvest_ohibc, file.path(dir_goal, 'data_explore/dfo_aq_harvest_ohibc.csv'))
DT::datatable(aq_harvest_ohibc)aq_plot_df <- aq_harvest_ohibc %>%
left_join(get_rgn_names(), by = 'rgn_id') %>%
mutate(species = str_replace(species, '_', ' '),
species = str_replace(species, '_[a-z_]+', ''),
species = ifelse(str_detect(species, 'sable'), 'sablefish', species))
aq_sum_df <- aq_plot_df %>%
group_by(year, rgn_name, aq_type) %>%
summarize(harvest = sum(harvest)) %>%
ungroup()
ggplot(aq_plot_df,
aes(x = year, y = harvest, color = species)) +
ggtheme_plot() +
geom_line(data = aq_sum_df %>% filter(aq_type == 'finfish'),
aes(x = year, y = harvest),
size = 1.5, linetype = 'dashed', color = 'grey30') +
geom_line(data = aq_sum_df %>% filter(aq_type == 'shellfish'),
aes(x = year, y = harvest),
size = 1.5, linetype = 'dotted', color = 'grey30') +
geom_point() +
geom_line(aes(group = species)) +
facet_wrap( ~ rgn_name, scale = 'free_y')Layers to ~/github/ohibc/prep/mar/v2017/output:
mapp_aq_types <- data.frame(rgn_id = c(1:4),
aq_type = c('shellfish', 'shellfish', 'finfish', 'shellfish'))
mapp_aq_areas <- read_csv(file.path(dir_goal, 'data_explore/mapp_aq_smz_areas.csv')) %>%
group_by(rgn_id, desc) %>%
summarize(area_km2 = sum(area_km2)) %>%
ungroup() %>%
left_join(mapp_aq_types, by = 'rgn_id') %>%
select(rgn_id, aq_type, area_km2) %>%
mutate(source = 'mapp')
dfo_aq_areas <- read_csv(file.path(dir_goal, 'data_explore/dfo_aq_tenure_areas.csv')) %>%
mutate(aq_type = ifelse(str_detect(filename, 'Finfish'), 'finfish', 'shellfish')) %>%
group_by(rgn_id, aq_type) %>%
summarize(area_km2 = sum(area_km2)) %>%
ungroup() %>%
distinct() %>%
mutate(source = 'dfo')
aq_areas <- bind_rows(mapp_aq_areas, dfo_aq_areas) %>%
group_by(rgn_id, aq_type) %>%
mutate(a_tot_km2 = sum(area_km2)) %>%
ungroup()
write_csv(aq_areas, file.path(dir_goal, 'output/aq_areas.csv'))
DT::datatable(aq_areas)To compare potential vs harvest, we need to convert bivalve units to metric tonnes. Some figures:
Averaging these gives about 27.5g per piece
aq_mass_per_pc <- 0.0275e-3 ### tonnes per piece
pot_b_unclipped <- read_csv(file.path(dir_goal, 'data_explore/prod_pot_b_range_unclipped.csv')) %>%
filter(prob_lvl == .50) %>%
select(rgn_id, median_prod_uncl = quants)
pot_b <- read_csv(file.path(dir_goal, 'data_explore/prod_pot_b_range.csv')) %>%
filter(prob_lvl == .50) %>%
select(rgn_id, median_prod = quants) %>%
left_join(pot_b_unclipped, by = 'rgn_id') %>%
mutate(median_prod_tonnes = ifelse(is.na(median_prod), median_prod_uncl, median_prod),
median_prod_tonnes = median_prod_tonnes * aq_mass_per_pc,
aq_type = 'shellfish',
units = 'tonnes') %>%
select(rgn_id, median_prod = median_prod_tonnes, aq_type, units)
pot_f <- read_csv(file.path(dir_goal, 'data_explore/prod_pot_f_range.csv')) %>%
filter(prob_lvl == 0.50) %>%
select(rgn_id, median_prod = quants) %>%
mutate(aq_type = 'finfish',
units = 'tonnes')
pot_aq <- bind_rows(pot_b, pot_f)
write_csv(pot_aq, file.path(dir_goal, 'output/aq_potential.csv'))
DT::datatable(pot_aq)harvest_df <- read_csv(file.path(dir_goal, 'data_explore/dfo_aq_harvest_ohibc.csv')) %>%
group_by(year, rgn_id, aq_type) %>%
summarize(harvest_tonnes = sum(harvest)) %>%
ungroup()
write_csv(harvest_df, file.path(dir_goal, 'output/aq_harvest.csv'))
DT::datatable(harvest_df)Reference points pictured are based only on DFO tenures; the MaPP reference areas are far larger and would create far more ambitious reference points for shellfish.
aq_area <- read_csv(file.path(dir_goal, 'output', 'aq_areas.csv'))
aq_potential <- read_csv(file.path(dir_goal, 'output', 'aq_potential.csv'))
aq_harvest <- read_csv(file.path(dir_goal, 'output', 'aq_harvest.csv'))
aq_all <- aq_area %>%
left_join(aq_potential, by = c('rgn_id', 'aq_type')) %>%
left_join(aq_harvest, by = c('rgn_id', 'aq_type')) %>%
mutate(ref_pt = area_km2 * median_prod) %>%
left_join(get_rgn_names(), by = 'rgn_id')
aq_f_df <- aq_all %>%
filter(aq_type == 'finfish') %>%
filter(source == 'dfo') %>%
select(year, rgn_name, harvest_tonnes, ref_pt)
aq_f_plot <- ggplot(aq_f_df, aes(x = year, y = harvest_tonnes)) +
ggtheme_plot() +
geom_line(aes(group = rgn_name)) +
geom_hline(aes(yintercept = ref_pt), color = 'red', alpha = .5) +
labs(title = 'Finfish harvest',
color = 'Ref pt source') +
facet_wrap( ~ rgn_name, scales = 'free_y') +
ylim(0, NA)
print(aq_f_plot)aq_b_df <- aq_all %>%
filter(aq_type == 'shellfish') %>%
filter(source == 'dfo') %>%
select(year, rgn_name, harvest_tonnes, ref_pt)
aq_b_plot <- ggplot(aq_b_df, aes(x = year, y = harvest_tonnes)) +
ggtheme_plot() +
geom_line(aes(group = rgn_name)) +
geom_hline(aes(yintercept = ref_pt), color = 'red', alpha = .5) +
labs(title = 'Shellfish harvest',
color = 'Ref pt source') +
facet_wrap( ~ rgn_name, scales = 'free_y') +
ylim(0, NA)
print(aq_b_plot)# prov_wrapup()